KRI <- read.csv("KRI.csv")
KRI$BLVC2F <- as.factor(KRI$BLVC2)
KRI$BHVC2F <- as.factor(KRI$BHVC2)
KRI$RESI <- as.factor(KRI$HOME)

library(lmtest)
library(sandwich)
library(modelsummary)
library(MASS)
library(ggplot2)

## Table 1
# Beliefs
mean(KRI$BLV1[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$BLV1[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$BLV1[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$BLV1[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$BLV1[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$BLV1[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$BLV1, na.rm=T)
sd(KRI$BLV1, na.rm=T)

mean(KRI$BLV2[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$BLV2[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$BLV2[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$BLV2[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$BLV2[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$BLV2[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$BLV2, na.rm=T)
sd(KRI$BLV2, na.rm=T)

mean(KRI$BLVC[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$BLVC[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$BLVC[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$BLVC[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$BLVC[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$BLVC[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$BLVC, na.rm=T)
sd(KRI$BLVC, na.rm=T)

t.test(KRI$BLV1[KRI$Atti=="Counter"], KRI$BLV2[KRI$Atti=="Counter"], na.rm=T)
t.test(KRI$BLV1[KRI$Atti=="Pro"], KRI$BLV2[KRI$Atti=="Pro"], na.rm=T)
t.test(KRI$BLV1[KRI$Atti=="Neutral"], KRI$BLV2[KRI$Atti=="Neutral"], na.rm=T)
t.test(KRI$BLV1, KRI$BLV2, na.rm=T)

# Intention
mean(KRI$SHARE1[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$SHARE1[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$SHARE1[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$SHARE1[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$SHARE1[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$SHARE1[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$SHARE1, na.rm=T)
sd(KRI$SHARE1, na.rm=T)

mean(KRI$SHARE2[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$SHARE2[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$SHARE2[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$SHARE2[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$SHARE2[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$SHARE2[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$SHARE2, na.rm=T)
sd(KRI$SHARE2, na.rm=T)

mean(KRI$SIC[KRI$Atti=="Counter"], na.rm=T)
sd(KRI$SIC[KRI$Atti=="Counter"], na.rm=T)

mean(KRI$SIC[KRI$Atti=="Pro"], na.rm=T)
sd(KRI$SIC[KRI$Atti=="Pro"], na.rm=T)

mean(KRI$SIC[KRI$Atti=="Neutral"], na.rm=T)
sd(KRI$SIC[KRI$Atti=="Neutral"], na.rm=T)

mean(KRI$SIC, na.rm=T)
sd(KRI$SIC, na.rm=T)

t.test(KRI$SHARE1[KRI$Atti=="Counter"], KRI$SHARE2[KRI$Atti=="Counter"], na.rm=T)
t.test(KRI$SHARE1[KRI$Atti=="Pro"], KRI$SHARE2[KRI$Atti=="Pro"], na.rm=T)
t.test(KRI$SHARE1[KRI$Atti=="Neutral"], KRI$SHARE2[KRI$Atti=="Neutral"], na.rm=T)
t.test(KRI$SHARE1, KRI$SHARE2, na.rm=T)


## Manipulation Check
t.test(KRI$BLV1[KRI$Atti=="Pro"], KRI$BLV1[KRI$Atti=="Counter"], na.rm=T)
t.test(KRI$BLV1[KRI$Atti=="Pro"], KRI$BLV1[KRI$Atti=="Neutral"], na.rm=T)
t.test(KRI$BLV1[KRI$Atti=="Counter"], KRI$BLV1[KRI$Atti=="Neutral"], na.rm=T)

#Models 1 to 10
KRI$Atti <- factor(KRI$Atti, levels = c("Pro", "Counter", "Neutral"))
KRI$Atti <- relevel(KRI$Atti, ref = "Neutral")

# Model 1
m1_beta <- lm(scale(BLVC) ~ Atti + scale(BLV1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) +  scale(AGE) + RESI, data=KRI)
m1r_beta <- coeftest(m1_beta, vcov=vcovHC(m1_beta, type="HC1"))
modelsummary(m1_beta, vcov="HC1", stars=TRUE) 

# Model 2
m2_beta <- lm(scale(BLVC) ~ DPRO + scale(BLV1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(AGE) + RESI, data=KRI)
m2r_beta <- coeftest(m2_beta, vcov=vcovHC(m2_beta, type="HC1"))
modelsummary(m2_beta, vcov="HC1", stars=TRUE) 

# Model 3
m3_beta <- lm(scale(SIC) ~ Atti + scale(BLVC) + scale(SHARE1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI)+ scale(SFREQ) + scale(AGE) + RESI, data=KRI)
m3r_beta <- coeftest(m3_beta, vcov=vcovHC(m3_beta, type="HC1"))
modelsummary(m3_beta, vcov="HC1", stars=TRUE) 

# Model 4
m4_beta <- lm(scale(SIC) ~ scale(DPRO) + scale(BLVC) + scale(SHARE1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(SFREQ) + scale(AGE) + RESI, data=KRI)
m4r_beta <- coeftest(m4_beta, vcov=vcovHC(m4_beta, type="HC1"))
modelsummary(m4_beta, vcov="HC1", stars=TRUE) 

# Model 5
m5_beta <- lm(scale(BLVC2) ~ Atti + scale(ALLTRUB1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI)  + scale(AGE) + RESI, data=KRI)
m5r_beta <- coeftest(m5_beta, vcov=vcovHC(m5_beta, type="HC1"))
modelsummary(m5_beta, vcov="HC1", stars=TRUE) 

# Model 6
m6_beta <- polr(BLVC2F ~ Atti + scale(ALLTRUB1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(AGE) + RESI , data=KRI, Hess=TRUE)
m6r_beta <- coeftest(m6_beta) 
modelsummary(m6_beta, stars=TRUE)

# Model 7
m7_beta <- lm(scale(BHVC2) ~ Atti + scale(BLVC2) + scale(ALLSHAREB1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(SFREQ) + scale(AGE) + RESI, data=KRI)
m7r_beta <- coeftest(m7_beta, vcov=vcovHC(m7_beta, type="HC1"))
modelsummary(m7_beta, vcov="HC1", stars=TRUE) 

# Model 8
m8_beta <- polr(BHVC2F ~ Atti + scale(BLVC2)+ scale(ALLSHAREB1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(SFREQ) + scale(AGE) + RESI, data=KRI, Hess=TRUE)
m8r_beta <- coeftest(m8_beta) 
modelsummary(m8_beta, stars=TRUE)


# Model 9
m9_beta <- lm(scale(BLVC) ~ scale(PRO) + scale(COUNTER) + scale(BLV1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(AGE) + RESI, data=KRI)
m9r_beta <- coeftest(m9_beta, vcov=vcovHC(m9_beta, type="HC1"))
modelsummary(m9_beta, vcov="HC1", stars=TRUE) 

# Model 10
m10_beta <- lm(scale(SIC) ~ scale(PRO) + scale(COUNTER) + scale(BLVC) + scale(SHARE1) + scale(INT) + scale(EDU) + SEXF + scale(FAMI) + scale(SFREQ) + scale(AGE) + RESI, data=KRI)
m10r_beta <- coeftest(m10_beta, vcov=vcovHC(m10_beta, type="HC1"))
modelsummary(m10_beta, vcov="HC1", stars=TRUE) 

## Figures
# Figure 3
coef_df1 <- data.frame(Variables = rownames(m1r_beta),
                            Estimate = m1r_beta[,"Estimate"],
                            Lower = m1r_beta[,"Estimate"] - 1.96 * m1r_beta[, "Std. Error"],
                            Upper = m1r_beta[,"Estimate"] + 1.96 * m1r_beta[, "Std. Error"],
                            Lower1 = m1r_beta[,"Estimate"] - m1r_beta[, "Std. Error"],
                            Upper1 = m1r_beta[,"Estimate"] + m1r_beta[, "Std. Error"])

names <- c("Counter-attitudinal", "Pro-attitudinal")
names <- factor(names, levels=names)
ggplot(coef_df1[c(3,2),], aes(x = names, y = Estimate)) + 
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  geom_linerange(aes(ymin = Lower1, ymax = Upper1), color="black", linewidth=1.2)+
  labs(x=" ", y="Estimated Effect in Standard Deviations") +
  theme_classic() + theme(axis.text = element_text(size=11),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ggtitle("DV: Belief Correction") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  ylim(-0.25, 0.05) +
  coord_flip()

# Figure 4
coef_df2 <- data.frame(Variables = rownames(m3r_beta),
                       Estimate = m3r_beta[,"Estimate"],
                       Lower = m3r_beta[,"Estimate"] - 1.96 * m3r_beta[, "Std. Error"],
                       Upper = m3r_beta[,"Estimate"] + 1.96 * m3r_beta[, "Std. Error"],
                       Lower1 = m3r_beta[,"Estimate"] - m3r_beta[, "Std. Error"],
                       Upper1 = m3r_beta[,"Estimate"] + m3r_beta[, "Std. Error"])
names2 <- c("Belief Correction", "Counter-Attitudinal", "Pro-Attitudinal")
names2 <- factor(names2, levels=names2)
ggplot(coef_df2[c(4,3,2),], aes(x = names2, y = Estimate)) + 
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  geom_linerange(aes(ymin = Lower1, ymax = Upper1), color="black", linewidth=1.2)+
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  labs(x=" ", y="Estimated Effect in Standard Deviations") +
  theme_classic() + theme(axis.text = element_text(size=11),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ylim(-0.15, 0.205) +
  ggtitle("DV: Intention Correction") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  coord_flip()

## Figure 5 (Top Left)
coef_df3 <- data.frame(Variables = rownames(m5r_beta),
                            Estimate = m5r_beta[,"Estimate"],
                            Lower = m5r_beta[,"Estimate"] - 1.96 * m5r_beta[, "Std. Error"],
                            Upper = m5r_beta[,"Estimate"] + 1.96 * m5r_beta[, "Std. Error"])

names <- c("Counter-attitudinal", "Pro-attitudinal")
names <- factor(names, levels=names)
ggplot(coef_df3[c(3,2),], aes(x = names, y = Estimate)) + 
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  labs(x=" ", y="Est.Effect in Std.Dev.") +
  theme_classic() + theme(axis.text = element_text(size=13),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ylim(-0.26, 0.12) +
  ggtitle("DV: Belief Correction (OLS)") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  coord_flip()

## Figure 5 sd (Top Right)
coef_df4 <- data.frame(Variables = rownames(m6r_beta),
                            Estimate = m6r_beta[,"Estimate"],
                            Lower = m6r_beta[,"Estimate"] - 1.96 * m6r_beta[, "Std. Error"],
                            Upper = m6r_beta[,"Estimate"] + 1.96 * m6r_beta[, "Std. Error"])

names <- c("Counter-attitudinal", "Pro-attitudinal")
names <- factor(names, levels=names)
ggplot(coef_df4[c(2,1),], aes(x = names, y = Estimate)) + 
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  labs(x=" ", y="Std. Coef. Est.") +
  theme_classic() + theme(axis.text = element_text(size=13),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ylim(-0.75, 0.4) +
  ggtitle("DV: Belief Correction (Ordered Logit)") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  coord_flip()


## Figure 5 (Bottom Left) 
coef_df5 <- data.frame(Variables = rownames(m7r_beta),
                            Estimate = m7r_beta[,"Estimate"],
                            Lower = m7r_beta[,"Estimate"] - 1.96 * m7r_beta[, "Std. Error"],
                            Upper = m7r_beta[,"Estimate"] + 1.96 * m7r_beta[, "Std. Error"],
                            Lower1 = m7r_beta[,"Estimate"] - m7r_beta[, "Std. Error"],
                            Upper1 = m7r_beta[,"Estimate"] + m7r_beta[, "Std. Error"])

names2 <- c("Belief Correction", "Counter-attitudinal", "Pro-attitudinal")
names2 <- factor(names2, levels=names2)
ggplot(coef_df5[c(4,3,2),], aes(x = names2, y = Estimate)) + 
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  geom_linerange(aes(ymin = Lower1, ymax = Upper1), color="black", linewidth=1.2)+
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  labs(x=" ", y="Est.Effect in Std.Dev") +
  theme_classic() + theme(axis.text = element_text(size=13),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ylim(-0.26, 0.12) +
  ggtitle("DV: Intention Correction (OLS)") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  coord_flip()

## Figure 6 (Bottom Right) 
coef_df6 <- data.frame(Variables = rownames(m8r_beta),
                             Estimate = m8r_beta[,"Estimate"],
                             Lower = m8r_beta[,"Estimate"] - 1.96 * m8r_beta[, "Std. Error"],
                             Upper = m8r_beta[,"Estimate"] + 1.96 * m8r_beta[, "Std. Error"],
                             Lower1 = m8r_beta[,"Estimate"] - m8r_beta[, "Std. Error"],
                             Upper1 = m8r_beta[,"Estimate"] + m8r_beta[, "Std. Error"])

names2 <- c("Belief Correction", "Counter-attitudinal", "Pro-attitudinal")
names2 <- factor(names2, levels=names2)
ggplot(coef_df6[c(3,2,1),], aes(x = names2, y = Estimate)) + 
  geom_linerange(aes(ymin = Lower, ymax = Upper),color = "black") +
  geom_linerange(aes(ymin = Lower1, ymax = Upper1), color="black", linewidth=1.2)+
  geom_point(size = 2.5, color = "black") + geom_hline(yintercept = 0, lty="dotted") +
  labs(x=" ", y="Std. Coef.") +
  theme_classic() + theme(axis.text = element_text(size=13),
                          panel.background = element_rect(fill = NA)) +
  geom_text(
    mapping = aes(label = round(Estimate, digits = 3)),
    hjust = -0.2,
    vjust = -0.5, size = 15/.pt) +
  ylim(-0.75, 0.4) +
  ggtitle("DV: Intention Correction (Ordered Logit)") + 
  theme(plot.title = element_text(hjust = 0.5, size=25)) +
  coord_flip()

## Cross-tabulations
# Dependent Variables
table(KRI$ALLTRUB1, KRI$Atti)
table(KRI$ALLTRUB2, KRI$Atti)
table(KRI$BLVC2F, KRI$Atti)
table(KRI$ALLSHAREB1, KRI$Atti)
table(KRI$ALLSHAREB2, KRI$Atti)
table(KRI$BHVC2F, KRI$Atti)

# Control Variables
table(KRI$SEXF, KRI$Atti)
table(KRI$AGE, KRI$Atti)
table(KRI$EDU, KRI$Atti)
table(KRI$HOME, KRI$Atti)
table(KRI$FAMI, KRI$Atti)
table(KRI$INT, KRI$Atti)
table(KRI$SFREQ, KRI$Atti)
